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Chapter  1 


INTRODUCTION 

In  recent  years  numerous  finite  element  models  based  on  the 
degenerate  solid  shell  concept  [l]  have  been  proposed  for  the  analysis 
of  thin  shell  structures  [2-4],  A  distinct  advantage  of  the  degenerate 
solid  shell  concept  is  that  it  can  be  used  for  finite  element  modeling 
of  arbitrary  shell  geometries  without  resorting  to  a  specific  shell 
theory.  However,  when  used  in  an  assumed  displacement  formulation,  the 
performance  of  a  degenerate  solid  shell  element  deteriorates  rapidly  as 
thickness  decreases.  This  phenomenon  is  called  locking  and  results 
from  the  inability  of  an  element  to  represent  a  zero  inplane  and 
transverse  strain  state  without  disrupting  the  bending  behavior  [5]. 

One  popular  method  used  in  an  attempt  to  alleviate  locking  is 
reduced/selective  integration  [6-11],  However,  reduced/selective 
integration  has  had  limited  success.  In  some  elements  the  locking 
effect  is  not  completely  alleviated  [12].  More  commonly,  when  the 
order  of  integration  is  reduced  sufficiently  to  alleviate  locking, 
spurious  kinematic  or  zero  strain  energy  modes  are  introduced  to  the 
element  stiffness  matrix.  These  kinematic  modes  can  be  suppressed  by 
adding  a  stabilization  matrix  to  the  stiffness  matrix  calculated  using 
reduced  integration  [13,14], 

An  alternate  method  suggested  by  Lee  and  Pian  [5]  is  to  use  a 
formulation  with  assumed  independent  strain  based  on  the 
Hel 1 i nger-Rei ssner  principle.  A  nine-node  shell  element  based  on  this 
approach  performed  very  well  on  a  variety  of  thin  shell  test  problems 


1 


[15].  This  formulation  provides  a  rational  mathematical  basis  for  the 
reduced  integration  scheme  [5,  16,  17], 

Recently,  Rhiu  and  Lee  successfully  developed  nine-node  and 
sixteen-node  degenerate  solid  shell  elements  using  a  new  mixed 
formulation  [3,  4]  which  is  also  based  on  the  Hel 1 i nger-Rei ssner 
principle  with  independent  strain.  In  this  case,  the  assumed 
independent  strain  is  divided  into  a  higher  order  and  a  lower  order 
part.  The  new  mixed  formulation  provides  a  rational  basis  for 
introducing  a  stabilization  matrix  to  the  reduced  integration  displace¬ 
ment  model . 

Another  approach  to  the  finite  element  analysis  of  shell  structures 
that  also  allows  for  the  modeling  of  arbitrary  geometries  without 
invoking  a  specific  shell  theory  is  to  use  a  three-dimensional, 
eighteen-node,  solid  element.  In  fact,  a  solid  element  is  more 
convenient  than  a  degenerate  solid  shell  element  since  it  does  not  need 
rotational  angles  to  describe  the  kinematics  of  deformation. 

However,  when  used  in  a  conventional  assumed  displacement  finite 
element  formulation  this  eighteen-node  element  performs  poorly.  It 
experiences  the  same  mplane  and  transverse  shear  locking  that 
degenerate  solid  shell  elements  do.  In  addition,  there  is  a  normal 
strain  locking  due  to  the  inability  of  the  element  to  represent  the 
condition  of  zero  normal  strain.  The  use  of  reduced  integration 
results  in  spurious  kinematic  or  zero  strain  energy  modes  in  the  finite 
element  model.  Hence,  the  global  stiffness  matrix  may  be  kinematically 
unstable.  Although  many  kinematically  unstable  elements  can  be  used 
with  care,  they  are  of  little  value  for  general  purpose  use. 


Therefore,  in  this  study  a  new  mixed  formulation  with  assumed 
independent  strain  [12]  is  applied  to  the  eighteen-node  solid  element 
to  improve  element  performance.  After  a  brief  discussion  of  geometry 
and  kinematics  of  deformation,  a  new  mixed  formulation  for  solid 
elements  is  presented.  Next,  the  selection  of  assumed  independent 
strain  is  discussed  in  detail.  This  is  followed  by  a  brief  discussion 
on  the  modified  stress-strain  relation.  Finally,  the  performance  of 
the  element  is  tested  by  solving  several  thin  plate  and  shell  problems 
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I 

GEOMETRY  AND  KINEMATICS 

i 

I 

I 

Figure  1  shows  the  isoparametric  representation  of  a  three-  j 

dimensional,  eighteen-node  element.  There  are  three  displacement  degrees 

I 

of  freedom  at  each  node.  Shape  function  polynomials  are  quadratic  in 
parent  coordinates  £  and  n,  and  linear  in  r,.  The  linear  coordinate  c 
is  in  the  thickness  direction  of  a  thin  structure. 

In  addition  to  a  global  Cartesian  coordinate  system  with  I 

I 

coordinates  X,  Y  and  Z,  local  orthogonal  coordinate  systems  are  used  to 
incorporate  shell  behavior  into  the  finite  element  formulation  as  will 
be  shown  later.  In  particular,  local  coordinate  systems  are  defin*  d  at 
numerical  integration  points.  This  enables  the  use  of  the  stress- 
strain  relation  as  well  as  strain  components  defined  with  respect  to 
the  local  coordinate  system.  For  a  local  coordinate  system,  the  three 
axes  x,  y  and  z  are  parallel  to  local  orthogonal  unit  vectors  a^, 
and  a-^,  respectively.  Local  coordinate  vector  a^  is  chosen  to  be 
parallel  with  either  £  or  n.  These  coordinate  systems  are  chosen  so 
that  the  element  stiffness  matrix  is  invariant  for  a  given  element 
geometry.  The  technique  for  choosing  coordinate  systems  at  integration 
points  is  discussed  in  detail  later. 

Given  these  coordinate  systems,  the  position  vector  X  in  the  global 
coordinate  system  of  an  arbitrary  point  is 


where  is  the  three-dimensional  Lagrange  shape  function  at  node  i  and 
X.  are  the  values  of  X  at  node  i.  Similarly,  the  global  displacement 


vector  U  of  the  same  point  is 


U  =  V  =  ^  Ni(.,n,t)yi 


where  is  the  vector  of  nodal  displacements. 


-G 


In  matrix  form,  the  linearized  engineering  strain  vector  e  with 


respect  to  global  coordinates  is 


U ,  v  +  V ,  s 


V.7  +  W„ 


u,7  +  w,. 


where  the  comma  (,)  denotes  differentiation.  For  a  finite  element 


formulation,  the  relationship  between  e  and  the  r.odal  degrees  of 


freedom  vector  a  can  be  written  symbolically  as 


7°  =  bg  9p 


P  — r 

where  B  is  a  matrix  relating  7’  to  a  .  The  strain  vector  in  the  local 


coordinate  system  is  found  by  transforming  the  global  strain  vector  as 


fol lows: 


7  =  T  7°  = 


kVa*!' 


where 


L£xx  eyy  ezz  exy  eyz  ezx-^ 


B  =  T  B 


In  Eq.  (5)  T  is  a  6x6  transformation  matrix  and  the  superscript  T 


denotes  the  transpose. 
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FINITE  ELEMENT  FORMULATION 


Assuming  small  displacements  and  no  initial  strains  (such  as 


thermal  strains),  the  functional  for  the  Hel 1 inger-Reissner 


principle  is 


If  (£T  £  £  -  j  £T  C  e)dV  -  M 


e  o 


where 


~  ~  ^-EXX  eyy  £zz  exy  £yz  ezx-* 


is  the  independent  strain  vector,  WQ  is  the  applied  load  term,  V0  is 


the  volume  of  the  element,  C  is  a  6x6  matrix  of  elastic  coefficients 


and  c  is  the  local  strain  vector  derived  from  the  displacement  field 


given  by  Eq.  (5).  The  summation  sign  indicates  assembly  of  all 


el ements. 


As  proposed  in  reference  12,  the  independent  local  strain  vector  is 


divided  into  two  parts  as  follows: 


£  =  £L  +  ~H 


where  and  are  the  independent,  local  strain  vectors  with  lower 


order  and  higher  order  assumed  polynomial  terms  in  £,  n  and  c , 


respectively.  Substituting  Eq.  (8)  into  Eq.  (6)  and  expanding  yields 


the  following  functional: 


(9) 


"R  =  n/ej  ~  ~  dve  "  7  ~  ~L  dVe  +  /&!  £  £ 


/  £h  £  ~L  dVe  "  1  ^  ~H  ~  ~H  dVe^  '  wo 


In  the  present  eighteen-node  element  the  lower  order  independent 
strains  are  assumed  to  be  trilinear  in  5,  n  and  5.  They  are  expressed 
in  terms  of  displacement-dependent  strains  evaluated  using  a  2x2x2 
point  Gauss  quadrature  rule.  Specifically, 


kU  /V  w 

~L  =  J-  Ni^  »n*d£i  =  J  Ni(5,n,t)B(ci,ni,?i)3e  (10) 

where  Ni  is  the  trilinear  shape  function  such  that  Npl  at  sampling 
point  i  of  the  2x2x2  point  integration  rule,  is  di spl acement- 
dependent  strain  evaluated  at  sampling  point  i,  and  r^,  C-j  are  the 
values  for  5,  n  and  t;  at  sampling  point  i.  Symbolically,  the  relation 
between  eL  and  qe  is  rewritten  as 

£L  =  I(C,n,c)qe  (11) 


where 


B(C,n,t)  =  l  N.(c,n,c)  B(5i,ni,ci) 


The  higher  order  strain  vector  Is  assumed  to  have  higher  order 
terms  in  £,  n  and  c ,  and  is  written  as 


£h  =  £(C,n,c)  <j 


k 

* 

I 


Buna? 


TV  Viwm  VWJ(  ’  VV* ITV^TiTl^rXTT!  .VV.V.VT.V.VT'VA-  V  v  V  v  ■■.’  V-  v 


where  P  is  the  shape  function  matrix  of  higher  order  assumed  strains 
and  a  is  a  vector  containing  element  strain  parameters. 

Introducing  Eqs.  (5),  (11)  and  (12)  into  Eq.  (9)  yields  the 
following  expression: 


where 


"r  =  I  (7  9e  +  21  5  3e  _  7  2T  2  “  9e  3e) 


£L  =  /(&T  3  8  +  BT  C  B  -  IT  C  B)  dVe 
G  =  /PT  C  B  dVg  -  /PT  C  J  dVg 
H  =  /PT  C  P  dV 


and  0  is  the  element  load  vector. 


Taking  5itR=0  with  respect  to  a  for  each  element  yields  the 
following  relationships: 

S  Se  -  H  a  =  0 


and 


2  =  a'1  s  % 


av 

(13a) 

(13b) 

(13c) 


(14a) 


(14b) 


Equation  (14a)  or  (14b)  represents  the  compatibility  equation  for  each 


(15b) 


and  the  stabilization  matrix  Ks  is  given  by 
K$  =  GT  H_1  G 

After  assembly  over  all  elements,  Eq.  (15)  can  be  rewritten  as 

>R  -  ?aT  £  a  -  sT  9  (i6) 

where  K  is  the  global  stiffness  matrix,  £  is  the  global  displacement 
vector  and  Q  is  the  global  load  vector.  Setting  6ttr  =  0  with  respect 
to  £  gives 

j<  3  =  9  (17) 


which  can  be  solved  for  £.  Thus,  £g  is  known,  and  the  local  strain 
vector  is  determined  as  follows: 


e  =  e,  +  P  G  q 

/v  'wl  >v  /v  /v  aJ 


(18a) 


Stress  a  is  determined  by  using  the  stress-strain  relation  and  the 
strain  in  Eq.  (18a)  as  follows: 


a  =  C  e 


(18b) 


As  will  be  demonstrated  in  the  next  section  the  J<L  matrix  in  Eq. 

(13a)  has  spurious  kinematic  modes.  The  kinematic  modes  occur  because 

only  lower  order  terms  of  assumed  strain  are  used  in  constructing  K.  . 

When  higher  order  assumed  strain  Is  properly  chosen,  the  K<~  matrix 

plays  the  role  of  a  stabilization  matrix  and  the  element  stiffness 

matrix  will  be  kinematically  stable. 

Before  discussing  the  selection  of  higher  order  strain  a  few 

comments  on  numerically  integrating  the  expressions  for  ir_  in  Eqs. 

K 
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xi 


(13a-c)  are  in  order.  Since  e  is  quadratic  in  5  and  n,  and  linear  in 
C  for  an  element  with  regular  parallelepiped  geometry,  and  is 
assumed  to  be  trilinear,  then  Eq.  (13a)  is  integrated  exactly  with  a 
2x2x2  point  Gauss  numerical  integration  rule.  In  this  case  Eq.  (13a) 
reduces  to 

K.  =  /  BT  C  B  dV  (19) 

~L  ^  ~  ~  ^  e 

where  the  subscript  L  on  the  integral  sign  indicates  a  2x2x2  point 
integration  rule.  In  other  words,  the  matrix  in  Eq.  (13a)  is 
exactly  the  same  as  the  stiffness  matrix  for  the  assumed  displacement 
model  with  a  2x2x2  point  reduced  integration  rule.  Detailed 
discussions  of  the  equivalence  between  reduced  integration  and  mixed 
formulations  are  in  references  12,  16  and  17.  Assuming  £H  is  at  most 
quadratic  in  5  and  n,  and  linear  in  c,  Eq.  (13c)  and  the  first  integral 
in  Eq.  (13b)  require  a  3x3x2  point  integration  rule  for  exact 
integration  of  the  same  regularly  shaped  element.  The  second  integral 
in  Eq.  (13b)  is  integrated  exactly  with  a  2x2x2  point  rule.  Although 
these  rules  are  only  exact  for  an  element  with  rectangular  geometry, 
the  same  integration  rules  are  adopted  for  elements  with  distorted 
geometries. 


•f. *.  ••  . 


Chapter  4 


HIGHER  ORDER  ASSUMED  STRAIN 

The  critical  step  in  the  present  formulation  is  the  proper  choice 
of  assumed  strain.  Generally,  the  assumed  strain  should  be  as  simple 
as  possible  to  avoid  locking.  However,  an  exceedingly  simple  assumed 
strain  field  will  trigger  spurious  kinematic  modes  that  do  not  produce 
strain  [3],  These  kinematic  modes  are  of  two  types;  compatible  and 
incompatible.  Compatible  kinematic  modes  persist  even  when  two  or  more 
elements  are  assembled.  In  other  words,  compatible  kinematic  modes  are 
modes  that  will  show  up  in  the  global  stiffness  matrix.  On  the  other 
hand,  incompatible  kinematic  modes  are  suppressed  when  two  or  more 
elements  are  a  sembled.  Since  incompatible  kinematic  modes  have  no 
adverse  effect  on  the  global  stiffness  matrix,  assumed  strain  can  be 
chosen  to  suppress  only  compatible  kinematic  modes  in  the  element 
stiffness  matrix.  In  addition,  assumed  strain  terms  must  be  chosen 
carefully  in  order  to  avoid  reintroducing  the  locking  effect. 


4.1  Local  Coordinate  System 

The  typical  result  of  choosing  assumed  strain  as  simple  as  possible 
is  an  incomplete  set  of  shape  function  polynomials  in  P.  In  general, 
this  leads  to  an  element  stiffness  matrix  which  Is  not  Invariant. 
Although  invariance  is  not  always  important,  it  can  be  enforced  by 
assigning  a  specific  local  coordinate  system  for  a  given  element 


geometry  [18].  For  an  arbitrary  thin  element,  two  unit  vectors  v1  and 
y2  are  defined  at£=n  =  c=  0so  that  v,  is  parallel  to  £  and  y2  is 


parallel  to  n.  The  angle  e  between  and  Vg  is  given  by 

9  =  cos"1  (v1  •  v^)  (20) 

If  9  is  less  than  or  equal  to  90°,  then  the  a^  unit  vector  in  the  x 
direction  of  the  local  coordinate  system  is  chosen  parallel  with  £. 
Otherwise  a^  is  chosen  parallel  to  n.  Next  the  a^  vector  normal  to  the 
5-n  plane  is  determined.  Finally,  unit  vector  ^  is  found  by  taking 
the  cross  product  between  a^  and  a^.  Rased  on  the  value  of  9  for  the 
element,  vectors  a^ ,  a^  and  a^  can  computed  at  any  point  in  the 
element.  In  particular,  they  are  calculated  at  each  numerical 
integration  point.  This  definition  of  a^,  a^  and  a3  guarantees  a 
unique  set  of  local  coordinate  systems  for  a  given  element  geometry, 
and  thus  a  stiffness  matrix  independent  of  choice  of  global  coordinate 
system. 

4.2  Spurious  Kinematic  Modes 

As  discussed  in  Chapter  3,  the  K|  matrix  in  Eq.  (13a)  or  (19)  has 
kinematic  modes  which  can  be  suppressed  with  the  proper  choice  of 
higher  order  assumed  strain.  Thus,  when  P  is  chosen  properly,  the  K<- 
matrix  piays  the  role  of  a  stabilization  matrix.  In  fact,  for  the  new 
mixed  formulation  the  higher  order  assumed  strain  is  chosen  specifi¬ 
cally  to  suppress  kinematic  modes  [12].  Therefore,  the  kinematic  modes 
must  be  known  before  the  higher  order  strain  terms  are  selected. 

Kinematic  modes  for  the  K^  matrix  can  be  determined  analytically 
for  an  element  with  regular  parallelepiped  geometry.  For  a  cubic 
element  with  sides  along  x  =  ±1,  y  =  ±1  and  z  =  ±1  and  the  global  and 


local  coordinate  systems  coincident,  the  displacement  field  can  be 
written  as 


where 
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The  local  displacement-dependent  strain  vector  is  determined  by 
introducing  Eq.  (21)  to  Eq.  (3).  The  displacement  field  corresponding 
to  zero  strain  is  found  by  setting  tiie  local  strain  vector  equal  to 
zero  at  the  2x2x2  integration  points.  Specifically, 

J  (il//3,tl//3,±l//3)  =0  (22) 

Equation  (22)  represents  48  homogeneous  equations  with  54  unknown 
variables.  Solving  these  equations  gives  the  following  displacement 


field  for  zero  strain: 

=  a i  >  a3y  +  a1Qz  +a2x(l  -  3y2 )  +  a4(x2  +  y2  -  3x2y2 ) 

+  auzx(l  -  3y")  +  a13z(x2  +  y2  -  3x2y2  )  (23a) 

V  =  b,  -  a3x  +  b1Qz  -  a2y(l  -  3x2 )  +  b4(x2  +  y2  -  3x2y2 ) 

-  a11zy(l  -  3x2 )  +  b13z(x2  +  y2  -  3x2y2 )  (23b) 
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W  —  c ^  -  3^qX  b i qV  —  1/3s^2^  ~l/3b^3y 

+  c4(^2  +  y2  -  3x2y2)  +  c1Qz(l  -  3x2  -  3y2  +  9x2y2 )  (23c) 

The  underlined  terms  in  Eqs.  ^3a-c)  are  rigid  body  modes.  The  other 
terms  represent  spurious  kinematic  modes. 

4.3  Assumed  Strain 

The  strain  terms  corresponding  to  the  kinematic  modes  in  Eqs.  ( 23 a-c ) 
are 

exx =  a2^  ■  3y2)  +  a4^2x  ■ 6**2)  +  dllz(1  "  3/  ) 


+  a13z(2x  -  6xy2  )  (24a) 

eyy  =  -a2(l  -  3x2)  +  b4(2y  -  6x2y)  -  anz(l  -  3x2  ) 

+  t>13z(2y  -  6x2y)  (24b) 

7ZZ  *  c1Q(l  -  3x2  -  3y2  +  9x2y2)  (24c) 

exy  =  a4^2-v  ~  6x2^  +  b4^2x  ‘  6x/)  +  al3z(2>'  ‘  6x2-/) 

+  b13z(2x  -  6xy2 )  (24d) 

eyz  =  ‘  all>'^1  "  3x2 )  +  bi3^_  J  +  x"  ¥  /  ~  3xV) 

+  c4(2y  -  6x2y)  +  c1Qz(-  6y  +  18x2y)  (24e) 

ezx  =  aiix(1  “  3>’2)  +  a 1 3 (  *  y  +  x2  +  y2  "  3x2y2) 

+  c4(2x  -  6xy2)  +  c1Qz(-  6x  +  18xy2  )  (24f) 


Higher  order  assumed  strain  terms  are  chosen  by  examining  Eqs. 
(24a-f).  For  example,  the  kinematic  mode  represented  by 


_  2  2 

exx  =  a4(2x-6xy  )  is  suppressed  if  xy  is  added  to  the  higher  order 
assumed  strain  shape  functions  for  UXX)H-  Since  the  kinematic  mode 
represented  by  a^  can  also  be  suppressed  by  including  an  x  y  term  in 
(eXy)H»  there  must  be  a  logical  basis  for  deciding  between 
the  options.  The  general  rule  is  to  minimize  the  total  number  of  terms 
in  the  assumed  strain  that  contribute  to  locking.  For  a  three- 
dimensional  element  applied  to  a  thin  structure  this  means  that  the 
polynomial  terms  independent  of  z  (thickness)  should  be  kept  as  simple 
as  possible.  Kinematic  modes  corresponding  to  a^  and  a ^  need  not  be 
suppressed  since  they  are  incompatible  and  disappear  for  the  assembly 
of  two  or  more  elements. 

Even  with  the  general  rule  there  are  still  two  reasonable  alternate 
sets  of  higher  order  strain.  One  option  is  to  minimize  the  total 
number  of  terms  in  the  higher  order  strain  field  as  follows: 


(exx>H  =  alxy2  +  a2^z 

( eyy ) H  =  a3x2*  +  a4*2yz 

=  °5XV 

^  G  xy  ^  H  =  0 
(£u,)u  =  <*(;X2y 


(25) 


'yz'H 
(ezx)H  *  a?xy2 

In  Eq.  (25),  a^,  a^,  a^,  a^,  a^,  a^,  a7  are  unknown  strain  parameters. 
Note  that  only  one  term  is  needed  between  (e^2)H  and  (ezx)H  to  suppress 
the  kinematic  mode  represented  by  c^.  The  second  term  is  added  to 
avoid  introducing  a  directional  imbalance  to  the  element  stiffness 


matrix.  The  model  using  this  choice  of  assumed  strain  is  designated 
New  Mixed  Solid,  18-node  version,  A  (NMS18A). 

The  other  alternative  is  to  replace  the  x2y2  term  in  (e  )u  with 

22  rl 

2  2 

x  yz  and  xy  z  terms  in  Uyz)H  and  UZX)H,  respectively.  Once  again 
only  one  term  is  needed  to  suppress  a  kinematic  mode  (in  this  case 
c^),  but  the  second  term  is  added  to  avoid  introducing  a  directional 
imbalance.  In  this  case,  the  total  number  of  terms  is  increased  by  one 
to  eight,  but  the  new  terms  are  expected  to  have  a  proportionally 
smaller  effect  as  the  thickness  becomes  smaller.  This  proposal  is 
tested  numerically  in  the  next  section.  This  formulation  is  designated 
NMS18B. 

For  completeness,  all  modes  are  suppressed  by  adding  to  the  NMS18B 
version  y  and  y  z  terms  to  UXX)H  and  x  and  x  z  terms  to  (eyy)H 
to  suppress  the  a?  and  incompatible  kinematic  modes.  As  in  the 
other  formulations,  only  one  term  is  needed  to  suppress  each  mode,  but 
two  terms  are  used  to  avoid  introducing  a  directional  imbalance.  This 
formulation  is  designated  NMS18C. 

For  an  element  with  arbitrary  geometry,  £,  n  and  ?  are  used  in 


place  of  x,  y  and  z  in  Eq.  (25).  Therefore,  with  the  local  coordinate 
system  defined  previously,  the  higher  order  assumed  strains  for  the 
eighteen-node  solid  element  are  chosen  as  follows: 


ITKm9  v*jr*jrr^jr*jr«jr«jn*jr*jrr' 


and  for  NMS18C 
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For  Eqs.  (27a),  (28a)  and  (29a),  f,  g,  h  and  i  are  chosen  as  follows: 
(1)  if  x  or  a^  is  parallel  to  c  (0  <  90°) 


til 


if  x  or  dj  is  parallel  to  n 
f  =  C2n 


(0  >  90°) 


9  =  Cn 


r  •  i 


s  =  n 


4.4  Modified  Stress-Strain  Relation 

The  stress-strai n  relation  for  an  isotropic,  three-dimensional 
solid  can  be  modified  to  incorporate  thin  plate  and  shell  behavior  by 
assuming  that  first  the  effect  of  a  on  e  and  e  is  small  and 

Z  Z  XX  Jr  Jr 

second  the  effect  of  a  and  a  on  c  is  small.  Then 

xx  yy  zz 
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xx  e  '  xx  yy; 
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(30a) 


(30b) 


(30c) 


where  E  is  Young's  Modulus  and  v  is  Poisson's  ratio.  The  relations 


between  shear  strains  and  stresses  remain  unchanged  as  follows: 

£xy  =  £  axy 


eyz  =  S  ayz 


C  >  V  0 

zx  G  zx 


where  G  is  the  shear  modulus. 


( 30d ) 


(30e) 


( 30f ) 


Inverting  Eqs.  (30a-f)  gives  the  modified  stress-strai n  relation  as 


fol lows: 


a  *  C_  e 
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Chapter  5 


NUMERICAL  TESTS 

Several  plate  and  shell  benchmark  problems  were  solved  to  test  the 
performance  of  the  eighteen-node  element  based  on  the  new  mixed 
formulation.  For  each  problem,  the  results  using  the  NMS18A  and  NMS18B 
formulation  are  compared  to  the  results  for  the  eighteen-node  assumed 
displacement  element  (designated  DISP18).  Results  for  the  NMS18C  model 
are  included  for  only  the  ring  problem  which  demonstrates  that  this 
model  locks  for  increasing  radius  to  thickness  ratios.  Results  using 
the  three-dimensional  and  modified  stress-strain  relations  are  compared 
for  plate  and  ring  displacement  problems.  In  addition,  analytical  or 
other  independent  solutions  are  presented  if  available.  All 
computations  were  done  in  double  precision  Fortran  on  the  UNIVAC 
1100/92  machine  at  the  University  of  Maryland. 

5.1  Simply  Supported  and  Clamped  Square  Plates 

A  square  flat  plate  under  uniform  pressure  is  a  good  starting  point 
for  testing  any  finite  element  that  is  to  be  used  for  shell  analysis. 
Effects  of  element  geometry,  boundary  conditions  and  length  to 
thickness  (L/t)  ratios  can  easily  be  compared  to  analytical  solutions. 
In  addition,  the  performance  of  the  £  and  matrices  is  easily  tested. 
The  deflections  of  a  square  plate  with  simply  supported  and  clamped 
edges  were  calculated  using  uniform  and  distorted  meshes.  Due  to 
symmetry  only  one  quarter  of  the  plate  is  modeled.  Typical  elastic 
properties  for  aluminum,  E  =  107  psi  and  v  =  0.3,  are  used. 


Evenly  divided  2x2,  3x3  and  4x4  uniform  meshes  were  used  for  both 
cases.  Distorted  2x2  and  4x4  meshes  were  adequate  for  the  simply 
supported  plate,  but  3x3  and  6x6  distorted  meshes  were  also  used  for 
the  clamped  plate  to  check  convergence.  The  6x6  distorted  mesh  is 
formed  by  bisecting  each  element  edge  in  the  3x3  distorted  mesh. 

Figures  2a  and  2b  depict  the  four  distorted  meshes.  Results  for  L/t 
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ratios  of  10  ,  10  and  10  are  presented. 

Normalized  maximum  deflections  for  uniform  meshes  are  in  Tables  1 
and  2;  distorted  mesh  results  are  in  Tables  3  and  4.  In  each  case  the 
maximum  deflection  at  the  centroid  of  the  plate  is  normalized  with 
respect  to  the  analytical  solution  obtained  from  thin  plate  theory 
[20,21].  The  NMS18A  and  NMS18B  formulations  give  numerical  results 
that  are  very  close  to  the  analytical  solution  over  a  wide  range  of  L/t 
ratios  when  uniform  meshes  and  the  modified  stress-strain  relation  are 
used.  NMS18C  results  are  not  reported,  since  for  plate  bending  NMS18C 
is  essentially  the  same  as  NMS18B.  This  is  because  the  two  models 
differ  only  in  inplane  strain  terms  and  there  is  no  inplane  locking  for 
the  plate  problem.  When  the  regular  stress-strain  relation  is  used, 
results  range  from  17  to  20  percent  below  the  analytical  solution. 

This  is  due  to  normal  strain  locking,  the  inability  of  the  element 
to  represent  accurately  the  condition  of  zero  normal  strain.  Results 
for  the  simply  supported  plate  with  distorted  meshes  and  the  modified 
stress-strain  relation  are  very  good.  The  clamped  plate  is  more 
sensitive  to  distorted  meshes  than  the  simply  supported  plate. 

However,  when  the  distorted  mesh  size  is  increased  to  6x6,  results  for 
the  clamped  plate  with  the  modified  stress-strain  relation  show  good 


agreement  with  the  analytical  solution.  Even  when  L/t  =  10  ,  which  is 
beyond  the  practical  range,  the  solution  differs  by  less  than  2.5%  from 
the  analytical  solution.  The  new  mixed  formulation  performed  better 
than  the  assumed  displacment  element  in  every  case.  The  assumed 
displacement  formulation  element  is  extremely  sensitive  to  increasing 
L/t  ratio  and  distorted  element  geometry.  Finally,  it  should  be  noted 
that  for  the  plate  problems,  the  NMS18A  and  NMS18B  elements  give 
virtually  identical  answers. 

Table  5  lists  the  nondimensional  bending  moments  for  a  4x4  uniform 
regular  mesh  evaluated  at  the  integration  point  nearest  to  the  centroid 
of  the  plate.  The  computed  values  are  normalized  with  respect  to  the 
analytical  solution  at  the  plate  centroid  using  thin  plate  theory 
[19,20,21].  Although  the  sampling  point  is  not  exactly  at  the 
centroid,  the  results  are  insensitive  to  changing  L/t  ratio.  This 
indicates  that  for  the  plate  case  the  new  mixed  formulation  elements 
are  reliable  for  stress  analysis. 

5.2  Pinched  Circular  Ring 

As  shown  in  Figure  3,  a  circular  ring  is  pinched  by  a  concentrated 
load  P  at  opposite  points  on  the  ring.  Due  to  symmetry  only  one 
quarter  of  the  ring  was  modeled  with  meshes  of  4,  8  and  16  equally 
divided  elements.  Material  and  geometric  properties  are  E  *  IQ7  psi, 
v  =  0.3  and  R/t  =  60,  100,  500. 

Table  6  lists  the  nondimensional  displacement  at  the  load  point 
normalized  to  the  analytical  solution  [5]  for  the  NMS18A,  NMS18B  and 
NMS18C  versions  using  the  regular  and  modified  stress-strai n 
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relationships.  No  results  for  the  assumed  displacement  model  are 
shown.  This  is  because  the  global  stiffness  matrices  resulting  from 
the  assumed  displacement  formulation  were  ill-conditioned  and  could  not 
be  solved  using  a  Cholesky  decomposition  scheme.  The  NMS18A  and  NMS18B 
formulations  give  accurate  results  for  all  mesh  sizes  when  the  modified 
stress-strain  relation  used.  The  results  for  the  NMS18A  and  NMS18B 
models  converge  to  a  value  about  9.5^  below  the  analytical  solution 
when  the  regular  stress-strain  relation  is  used.  NMS18C  exhibits 
locking  as  the  R./t  ratio  increases.  In  addition,  it  performs  more 
poorly  than  NMS18A  or  NMS18B  for  the  4  and  8  element  meshes. 

5.3  Pinched  Cylindrical  Shell 
5.3.1  Diaphragmed  Ends 

A  pinched  cylindrical  shell  with  diaphragmed  ends  is  a  good  deep 
shell  test  problem  since  an  analytical  solution  is  available  for 
comparison.  The  cylinder  is  pinched  by  a  concentrated  load  P  at  two 
opposing  points  on  the  circle  at  the  midsection.  Due  to  symmetry  and 
loading  only  one  octant  of  the  cylinder  was  modeled  (Figure  4). 

Uniform  and  refined  4x6,  5x7,  6x8  and  7x9  meshes  were  used.  The  4x6F 
(refined)  mesh  is  formed  by  dividing  the  elements  along  lines  BC  and  CD 
of  a  3x5  (uniform)  mesh.  The  3x5  and  4x6F  meshes  are  shown  in  Figure 
5.  The  meshes  illustrated  in  Figure  5  are  on  the  stretched  plane  of 
octant  ABCD  of  the  cylinder.  Refined  meshes  are  more  effective  in 
modeling  the  steep  gradients  of  deflections  and  stresses  near  the  con- 
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concentrated  load.  Material  and  geometric  properties  are  E  =  1.05x10 


Tables  7  and  8  list  the  nondimensional  displacement  Wc  =  -WcEt/P  at 
the  load  point  for  the  DISP18,  NMS18A  and  NMS18B  models  with  the 
modified  stress-strain  relation.  Results  using  the  regular  stress- 
strain  relation  are  not  reported  since  the  plate  and  ring  problems 
demonstrate  that  when  it  is  used  locking  occurs.  Since  the  NMS18C 

formulation  performed  poorly  on  the  ring  problem  no  further  results  for 

it  are  reported. 

Analytical  solutions  were  calculated  using  100  terms  in  each 

direction  for  the  double  Fourier  series  expression  given  by  Flugge  [22] 

and  reported  in  Reference  19.  For  comparison,  the  results  for  a  new 
mixed  formulation  degenerate  solid  shell  element  are  reported.  The 
element  is  the  SHEL9N  element  that  was  first  investigated  by  Rhiu  [19] 
and  examined  in  further  detail  in  References  12  and  18.  The  SHEL9N 
element  has  nine  nodes  with  three  displacement  and  two  rotational 
degrees  of  freedom  at  each  node.  Therefore,  SHEL9N  has  45  degrees  of 
freedom  versus  54  for  the  present  eighteen-node  solid  elements. 

Once  again,  the  assumed  displacement  model  performs  poorly, 
especially  for  the  larger  R/t  ratios.  For  R/t  =  100,  both  the  NMS18A 
and  NMS18B  results  are  close  to  the  analytical  solution  and  compare 
favorably  with  SHEL9N.  When  the  R/t  =  500  case  is  examined,  it  becomes 
clearer  that  the  NMS18B  formulation  is  superior  to  the  NMS18A 
formulation.  In  fact,  the  NMS18B  results  are  virtually  identical  to 
the  SHEL9N  answers. 

5.3.2  Clamped  Ends 

For  this  test  case  the  cylinder  of  the  previous  example  is  used. 


but  the  ends  are  clamped  instead  of  diaphragmed.  Results  for  two 

degenerate  solid  shell  elements,  one  a  triangular  element  [2]  and  the 

other  a  nine-node  element  based  on  a  conventional  mixed  formulation 

[15],  are  presented  for  comparison.  Nondimensional  deflections  at  the 

load  point  C  are  listed  in  Tables  9  and  10  for  R/t  =  100  and  500, 

respectively.  As  before,  the  assumed  displacement  element  performs 

poorly.  And  once  again,  while  answers  for  NMS18A  and  NMS18B  are  very 

close  at  R/t  =  100,  NMS18B  is  much  better  at  R/t  =  500.  In  addition. 

Figures  6  and  7  show  inplane  force  N  and  moment  M  per  unit  length 

x  y 

along  line  BC  for  a  9x7  mesh  with  NMS18B  elements.  The  results  in 
Figures  6  and  7  show  excellent  agreement  with  reference  4.  The  finite 
element  used  in  reference  4  is  an  accurate,  cubic  degenerate  solid 
shell  element.  Results  for  this  element  are  not  shown  to  avoid 
cluttering. 

5.4  Hemisphere  under  Alternating  Point  Loads 

A  hemishpere  under  directionally  alternating  point  loads  at  the 
free  edge  is  shown  in  Figure  3.  Due  to  symmetry  only  one  quarter  of 
the  hemisphere  was  modeled.  In  order  to  use  only  eighteen-node 
elements,  a  0.1°  cutout  was  made  at  the  pole  and  inplane  displacements 
along  the  resulting  edge  were  constrained  to  zero.  Both  uniform  and 
refined  4x5,  5x6,  6x7  and  7x8  meshes  were  used.  For  the  uniform  meshes 
each  element  subtends  equal  angles  in  the  longitude  (<j>)  and  colatitude 
(9)  directions.  For  a  refined  mesh  the  row  of  elements  nearest  the 
pole  are  divided  equally  in  the  colatitude  direction.  In  this  way  a 
4x5F  (refined)  mesh  is  formed  from  a  4x4  uniform  mesh.  Figure  8  shows 


a  4x4  uniform  mesh  divided  to  form  a  4x5  mesh.  For  convenience  the 
cutout  is  made  after  the  mesh  is  formed.  Material  and  geometric 
constants  are  E  =  107  psi,  v  =  0.3,  R  =  10  in.,  P  =  2  and  R/t  =  250, 
500. 

Nondimensional  displacement  TTA  =  DWA/PR2  at  point  A  is  reported  in 
Tables  11  and  12  for  R/t  =  250  and  500,  respectively.  The  symbol  0 
represents  bending  rigidity.  As  expected,  the  assumed  displacement 
model  performs  poorly.  For  R/t  -  250,  NMS18B  converges  more  quickly 
and  both  NMS18A  and  NMS18B  give  answers  close  to  Morley's  analytical 
solution  [23].  Morley  did  not  consider  the  R/t  =  500  case,  but  a 
converged  finite  element  solution  using  a  sixteen-node  element  similar 
to  SHEL9N  is  presented  for  comparison.  Results  for  the  three 
formulations  are  very  close. 
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Chapter  6 


CONCLUSIONS 


The  numerical  tests  demonstrate  that  the  NMS13A  and  NMS18B  mixed 
formulation  elements  with  the  modified  stress-strain  relation  give 
reliable  solutions  for  thin  plate  and  shell  problems.  In  the  present 
formulation ,  the  assumed  strain  in  conjunction  with  the  modified 
stress-stra i n  relation  effectively  eliminates  inplane,  shear  and  normal 
strain  locking.  Overall,  the  NMS18B  version  performed  better  than  the 
NMS13A  version,  thereby  supporting  the  proposal  that  higher  order 
strain  terms  with  a  thickness  coordinate  are  less  likely  to  reintroduce 
locking  than  other  terms.  Both  the  NMS13A  and  NMS13B  formulations  are 
kinematically  unstable  at  the  element  level,  but  are  stable  when  two  or 
more  elements  are  assembled.  The  NMS18C  formulation,  although  kinema¬ 
tically  stable  at  the  element  level,  performed  poorly  for  curved 
geometries.  Thus,  a  judiciously  chosen  higher  order  strain  has 
successfully  suppressed  compatible  kinematic  modes  without 
reintroducing  the  locking  effect.  Finally,  the  kinematics  of 
deformation  are  more  easily  described  with  the  eighteen-node  solid 
element  than  a  degenerate  solid  shell  element.  Therefore,  to  make  the 
most  of  this  advantage,  the  present  formulation  should  be  extended  to 
geometrically  and  materially  nonlinear  problems. 
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By 


Figure  2a  Distorted  2x2  and  4x4  mexhes  for  one  quarter 
of  a  square  plate 
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R/t=  60,100,500 
E  =  I07  psi 
v-O.Z 
P=  1.0 


Figure  3  Pinched  circular  ring 


Normalised  maximum  deflection  of  a  simply  supported  square  plate 
under  uniform  pressure  (regular  meshes) 


L/t 

Type 

2x2 

3x3 

4x 

A 

'1 

C 

Cm 

C 

C 

C 

C 

~m 

/V 

~m 

~m 

DISP13 

.7815 

.9573 

.8024 

.9830 

.8095 

.9918 

102 

NMS18A 

.8145 

1.0031 

.8170 

1.0033 

.8183 

NMS18B 

.8145 

1.0031 

.8169 

1.0033 

.8183 

DISP13 

.7802 

.9557 

.8008 

.9810 

.8077 

.9895 

103 

NMS18A 

.8135 

1.0018 

.8153 

1.0012 

.8159 

1.0008 

NMS18B 

.8135 

1.0018 

.8153 

1.0012 

.8159 

1.0008 

DISP18 

.7850 

.9585 

.8034 

.9846 

.8099 

.9919 

io4 

NMS18A 

.8135 

1.0013 

.8143 

1.0003 

.8147 

.9986 

NMS18B 

.8134 

1.0013 

.8143 

1.0006 

.8148 

.9992 

DISP13 
102  NMS18A 
NMS18B 


0ISP18 


D I SP 18 
NMS13A 


S 


2x2 

3x3 

C 

C 

C 

C 

/V 

~m 

.6643 

.7877 

.7427 

.9084 

.8271 

1.0129 

.8200 

1.0041 

.8271 

1.0129 

.8200 

1.0041 

.6373 

.7807 

.7362 

.9018 

.8254 

1.0111 

.8183 

1.0024 

.8254 

1.0111 

.8183 

1.0024 

.6383 

.7812 

.7368 

.9027 

.8255 

1.0111 

.8179 

1.0022 

.8255 

1.0112 

.8180 

1.0023 

.7712  .9448 
.8171  1.0009 
.8171  1.0009 


8  .945 


8168 

1.0004 

8168 

1.0005 

Table  3 


Normalized  maximum  deflection  of  a  simply  supported  square  plate 
under  uniform  pressure  (distorted  meshes) 


L/t 

Type 

2x2 

4x4 

0ISP18 

.9515 

.9930 

102 

NMS13A 

1.0260 

1.0129 

NMS18B 

1.0260 

1.0129 

DISP18 

.8757 

.9670 

103 

NMS18A 

1.0236 

1.0089 

NMS18B 

1.0236 

1.0089 

D I SP 13 

.1690 

.5119 

104 

NMS18A 

1.0234 

1.0059 

NMS18B 

1.0235 

1.0061 

!  L/t 

Type 

2x2 

3x3 

4x4 

6x6 

► 

: 

DISP18 

.5666 

.8394 

.9046 

.9727 

1 

[  102 

NMS18A 

.9432 

.9914 

1.0062 

1.0003 

m 

1 

NMS18B 

.9432 

.9914 

1.0062 

1.0003 

j  103 

DISP18 

NMS18A 

.0625 

.6117 

.1648 

.8978 

.6838 

.9551 

.9340 

.9935 

1 

I 

NMS18B 

.6117 

.8978 

.9551 

.9935 

1 

DISP18 

.0007 

.0020 

.0433 

.3815 

V, 

y 

104 

[ 

NMS18A 

.3623 

.8414 

.8978 

.9771 

» 

} 

NMS18B 

.3623 

.8415 

.8977 

.9772 

r.*, 

« + 

/- 

1 

: 

i 

» 

> 

i 

i 

i 

* 

42 

* 

■x 

$ 

1 

Bending  moment  Mx  at  integration  point  nearest  to  the  centroid 
of  a  square  plate  normalized  to  the  analytical  solution  at  the  centroid 

(uniform  4x4  mesh) 


Simply  Supported  Plate 

Clamped 

Plate 

L/t 

NMS18A 

NMS18B 

NMS1SA 

nmsibb 

102 

.9982 

.9982 

.9909 

.9909 

103 

.9961 

.9961 

.9907 

.9907 

io4 

.9938 

.9944 

.9901 

.9903 

Type 

4  elements 

8  elements 

c  o 

~  ~m 

C  C 

~  ~m 

NMS18A 

.8988 

1.0016 

.9157 

1.0043 

NMS18B 

.9145 

1.0029 

.9156 

1.0042 

NMS13C 

.6801 

.8940 

.8953 

.9895 

NMS13A 

.8895 

1.0009 

.9153 

1.0041 

NMS18B 

.9142 

1.0024 

.9154 

1.0040 

NMS1SC 

.5426 

.7648 

.8688 

.9648 

NMS18A 

.8570 

.9991 

.9146 

1.0035 

NMS13B 

.9134 

1.0008 

.9152 

1.0035 

NMS18C 

.2958 

.1785 

.7089 

.5542 

.9157  1.0044 
.9156  1.0044 
.9106  .9920 


.9154  1.00' 
.9154  1.00 
.8865  .88 


Nondimensional  deflection  at  point  C  (Wc=  -W  Et/P)  for  a 
pinched  cylinder  with  clamped  ends  (R/t  -  100) 


Type 

DISP18 

NMS18A 

NMS18B 

Mesh 

1 

Uniform 

Refined 

Uniform 

Refined 

Uniform 

Rof ined 

43.02 

74.75 

134.3 

135.7 

135.1 

135.8 

57.24 

88.18 

135.4 

136.3 

135.7 

136.2 

| 

70.66 

99.29 

136.0 

136.7 

136.0 

136.5 

82.29 

108.3 

136.4 

136.9 

136.3 

136.7 

Reference  2  137.01 


Reference  15  136.81 
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Table  10 


Nondimensional  deflection  at  point  C  (W  =  -WcEt/P)  for  a 


pinched  cylinder  with  clamped  ends  (R/t  =  500) 


DISP18 

NMS18A 

Uni  form 

Refined 

Uniform 

Refined 

49.37 

136.0 

784.2 

932.1 

72.59 

201.8 

859.4 

946.4 

101.5 

271.4 

906.0 

958.4 

135.3 

338.5 

933.1 

966.1 

Uni  form 


894.5 

928.6 
949.4 

960.6 


Refined 


963. 

964. 
967. 
971. 


Reference  2  963.93 


Reference  15  960.88 


W5S5 
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Table  11 

Nondimensional  deflection  at  point  A  (W^DW^PR2)  for  a  hemisphere 
under  alternating  point  loads  (R/t  =  250) 


D I SP 18 

NMS18A 

Regul ar 

Refined 

Regular 

Refined 

.0040 

.0030 

.1767 

.1799 

.0090 

.0069 

.1828 

.1845 

.0171 

.0135 

.1845 

.1855 

.0284 

.0233 

.1851 

.1857 

NMS18B 


Regular 

.1783 

.1834 

.1848 

.1853 


Ref i ned 
.1815 
.1351 
.1858 
.1858 


Nondimensional  deflection  at  point  A  (W^=DWA/PR2)  for  a  hemisphere 
under  alternating  point  loads  (R/t  =  500) 
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